Naval  Research  Laboratory 

Stennis  Space  Center,  MS  39529-5004 


NRL/FR/7440~03-10,053 


Wave  Bottom  Boundary  Layer  Models  for 
Smooth  and  Rough  Beds 


Jack  Puleo 

Mapping,  Charting,  and  Geodesy  Branch 
Marine  Geosciences  Division 

Olec  Mouraenko 

University  of  Florida 
Gcdnesville,  Florida  3261 1-6550 


September  17, 2003 


20031015  027 


Approved  for  public  release;  distribution  is  unlimited. 


15.  SUBJECT  TERMS 


wave  boiiom  turbulent  eddy  \'iscosiry 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

18.  NUMBER 

19a.  NAME  OF  RESPONSIBLE  PERSON 

OF  ABSTRACT 

OF  PAGES 

Jack  Puleo 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

UL 

22 

19b.  TELEPHONE  NUMBER  (Made  area 

Unclassified 

Unclassified 

Unclassified 

228-688-4328 

Standard  Form  298  (Rev.  8-98) 
l>rMcrib«l  by  ANSI  Std.  Z3ai6 

i 


CONTENTS 

INTRODUCTION . 1 

BOUNDARY  LAYER  THEORY  AND  TURBULENCE  MODELS . 2 

Turbulent  Eddy  Viscosity  Specification . 2 

One-Equation  Model . . . -3 

Two-Equation  k-z  Model . 4 

Two-Equation  Ar-co  Model . 5 

NUMERICAL  METHODS . -5 

Finite  Difference  Equations . 5 

Boundary  Conditions . 6 

Free-stream  Forcing . 2 

Boundary  Layer  Estimation . - . 7 

SOLUTION  PROCEDURE . 8 

MODEL  RESULTS . 9 

MODEL-DATA  COMPARISON... . . . . . 1 2 

DISCUSSION  AND  CONCLUSIONS . 17 

ACKNOWLEDGMENTS . 17 

REFERENCES . 18 

APPENDIX — Matlab  Codes . ; . . . 1 9 


iii 


WAVE  BOTTOM  BOUNDARY  LAYER  MODELS  FOR  SMOOTH  AND  ROUGH  BEDS 


INTRODUCTION 


The  goal  of  the  6.1  Foreshore  Sediment  Transport  Study  is  to  develop  predictive  capabilities  for  the 
hydrodynamics  and  sediment  transport  in  the  inner  surf  and  swash  zones  using  field  data  and  numerical 
models.  The  sediment  transport  aspect  of  the  study  often  stems  fi’om  some  prediction  of  the  fluid-supplied 
bed-shear  stress  resulting  from  fiiction  and  turbulent  fluctuations  near  the  bed  in  what  is  termed  the  wave 
bottom  boimdary  layer  (WBBL).  Hence,  knowledge  of  the  time-dependent  nature  of  fluid  flows  in  this 
region  will  lead  to  better  predictions  of  the  turbulent  motions  and  shear  stresses  and  ultimately  assist  in 
improving  sediment  transport  formulations  in  wave-dominated  environments. 


Much  of  the  work  on  WBBLs  has  focused  on  developing  various  forms  for  the  vertical  distribution  of 
the  turbulent  eddy  viscosity  (generically  referred  to  here  as  eddy  viscosity  models),  an  approximation 
assuming  that  turbulent  motions  diffuse  momentum  in  a  manner  similar  to  that  of  a  molecular  diftusive 
process  although  generally  at  a  much  faster  rate.  The  form  of  the  vertical  distribution  has  varied  from 
linear  to  parabolic  to  combined  forms  that  do  not  have  a  single  analytic  function  describing  them.  For 
instance,  Fredsoe  and  Deigaard  (1992)  discuss  several  forms  for  tiie  eddy  viscosity,  including  those  that 
can  be  ad^ted  to  independently  account  for  a  superimposed  current. 


A  more  detailed  approach  can  be  used  by  trying  to  calculate,  rather  than  model,  the  turbulent  motions 
themselves.  For  instance,  Spalart  (1988)  and  Moneris  and  Slinn  (submitted)  used  a  direct  numerical 
simulation  (DNS),  such  that  dl  scales  of  turbulent  motion  were  calculated,  rendering  a  description  of  the 
turbulent  eddy  viscosity  unnecessary.  While  DNS  simulations  are  probably  the  most  accurate  for 
determining  Ae  fluid  motions  in  the  WBBL,  their  extensive  computation  requirements  make  them 
unusable  for  rapid  estimations'  of  the  turbulent  fluid  motions  in  this  region. 


Rather  than  specifying  the  sh^e  of  the  turbulent  eddy  viscosity  or  calculating  the  fluid  motions  at  all 
scales,  a  middle  ground  approach  can  be  taken  where  the  turbulent  motions  are  modeled  leading  to  a  time 
and  space  varying  eddy  viscosity  (generically  referred  to  here  as  turbulence  closure  models).  One 
approach  is  referred  to  as  the  one-equation  or  k  (turbulent  kinetic  energy)  model.  Here  the  turbulent 
kinetic  energy  fluctuations  are  modeled  and  the  eddy  viscosity  is  estimated  by  specifying  some  length 
scale  of  the  turbulence  based  on  the  elevation  above  tte  bed  (e.g.,  Wilcox  2000).  A  second  approach  is  to 
use  a  two-equation  model  where  the  length  scale  is  no  longer  held  constant  and  itself  is  allowed  to  vary 
(e.g.,  Rodi  1980,  Justesen  1988,  and  Wilcox  1988). 

The  purpose  of  this  report  is  to  analyze  seven  WBBL  models  incorporating  various  methods  for 
estimating  the  turbulent  eddy  viscosity  that  can  be  used  to  rapidly  predict  the  fluid  flow  structure  and 
shear  stress  occurring  near  fte  bed.  The  models  will  be  compared  using  two  test  cases  and  against  a 
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laboratory  data  set.  The  report  includes  an  appendix  with  the 
description  of  how  to  run  the  simulations. 


computer  codes  used  (on  disk),  and  a  short 


BOUNDARY  LAYER  THEORY  AND  TURBULENCE  MODELS 


The  time-dependent  equatiori  of  turbulent  motion  for  an  incompressible 
ensemble  mean  oscillatory  velocity  component «  parallel  to  the  bed  is  given  by 


Newtonian  fluid  with 


an 


du 

aT 


p  dx  dz 


(1) 


where  p  is  the  fluid  density,  p  is  the  pressure,  v  is  the  kinematic 
ed(fy  viscosity,  x  is  the  horizontal  coordinate,  and  z  is  the  vertical 


viscosity  of  the  fluid,  o,  is  the  turbulent 
coordinate. . 


Outside  die  bound^  layer,  in  the  free  stream,  shear  stresses  are  no  longer  present  such  that  there  is 
balance  between  the  time-dependent  free-stream  velocity  and  the  horizontal  pr^sure  gradient  as 


dU  1  dp 

dt  pdx’  (2) 


where  C/is  the  free-stream  velocity. 


Substituting  Eq.  (2)  into  Eq.  (1)  and  rearranging  yields 


dt 


dz 


(u  +  y,) 


dz 


(3) 


diat  of^  ^  difference  between  the  depth-dependent  velocity  and 

reouire  O'’  specify  v,.  One  option  is  to  assume  or 

S  <■”  Many  other  opdons  exia,  several  of 


Torbulent  Eddy  Viscosity  Speciflcation 

Borrowing  from  open  channel  flow,  one  can  specify  the  shape  of  the  tuitulent  eddy  viscosity  a  priori 
^  some  function  d^nding  on  the  height  above  the  bed  and  the  friction  velocity  at  the  bed.  -^e  height 

turbulent  ^dy  size  is  restricted  as  the  bed  is  approached.  Friction  velocfly 
dependence  anses  to  incorporate  the  shear  at  the  bed  (and  hence  a  measOre  of  the  strength  of  the  flow). 

The  linear  form 


Ui  =  KU^Z, 


(4) 


3 


Wave  Bottom  Boundary  Layer  Models  for  Smooth  and  Rough  Beds 

where  k  is  the  von  Karman’s  constant  (=  0.4)  and  «•  is  the  friction  velocity  was  first  used  by,  for  example, 
Jonsson  (1966),  Kajiura  (1968),  and  Grant  and  Madsen  (1979),  allowing  for  an  analytic  solutioii  to 
Eq.(3). 

The  friction  velocity  is  given  as 


and  the  shear  stress  x  is,  in  turn,  defined  as 

T  =  p{v  +  o,)^.  (6) 

oz. 

Equation  (6)  is  often  seen  without  the  kinematic  viscosity  term,  however  we  include  it  since  the  flows 
will  typically  begin  from  rest  and  some  form  of  shear  near  the  bed  will  still  exist  on  the  molecular  level. 

The  linear  form  gives  the  expected  behavior  of  the  turbulent  eddy  viscosity  near  the  bed,  but  is 
physically  incorrect  near  the  surface.  The  turbulent  eddies  near  the  surface  are  predicted  to  be  too  large 
according  to  Eq.  (4),  since  it  is  expected  that  turbulent  shear  stresses  vanish  outside  the  boundary  layer 
(assuming  no  mean  currents  are  superimposed). 

To  relieve  some  of  this  difficulty,  Gelfenbaum  and  Smith  (1986)  modified  the  linear  turbulent  eddy 
viscosity  by  including  an  exponential  term  that  decays  with  increasing  elevation  above  the  bed.  Beach 
and  Sternberg  (1988)  further  modified  the  linear-exponential  form  to 

O,  ,  (7) 

where  D  is  the  domain  depth  in  their  study  of  sediment  suspension,  assuming  that  the  diffusion  of  mass 
and  momentum  were  equal. 

Finally,  Fredsoe  and  Deigaard  (1992)  give  a  parabolic  form  of  the  turbulent  eddy  viscosity  as 

Vf=KU*Z^-^^.  (8) 

It  may  expected  that  this  latter  form  would  be  more  correct  compared  to  Eqs.  (7)  and  (4)  since  it  has 
roughly  the  same  sh^e  near  the  bed  but  goes  to  zero  at  the  surface.  Equation  (8)  can  also  be  modified 
such  that  D  is  replaced  with  the  boundary  layer  height. 

One-Equation  Model 

Rather  than  specifying  the  form  of  the  eddy  viscosity,  it  can  be  estimated  based  on  the  turbulent 
kinetic  energy  in  the  flow  field.  The  turbulent  kinetic  energy  k  for  a  one-dimensional  problem  is  defined 
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where  uj  are  the  turbulent  fluctuations  from  the  mean  velocity.  The  transport  equation  for  k  is  derived  in, 
for  example,  Fredsoe  and  Deigaard  (1992),  and  given  as 


dt 


dz 


(u+u,) 


dz 


+  o, 


dud 

dz 


-C, 


(10) 


where  C;  is  a  constant  (=  0.09)  and  Z  is  a  turbulence  length  scale  given  by 


(11) 


In  Eq.  (10),  the  time  rate  of  change  of  turbulent  kinetic  energy  results  from  diffusion  (both  molecular 
and  turbulent),  production  (note  there  is  no  molecular  production  of  turbulent  kinetic  energy,  hence  v  is 
not  included  in  this  term),  and  dissipation.  Finally,  the  eddy  viscosity  results  from  the  Kolmogorov- 
Prandtl  relationship  between  eddy  viscosity,  a  mixing  length  and  the  turbulent  kinetic  energy  as 


(12) 


Two-Equation  k-£  Model 

Similar  to  the  one-equation  model,  the  turbulent  kinetic  energy  is  modeled  through  an  equation  similar 
to  Eq.  (10)  in  the  two-equation  model.  The  turbulence  length  scale,  however,  is  allowed  to  vary  in  time 
and  space  rather  than  remain  fixed  as  in  the  one-equation  model.  In  order  for  L  to  vary,  an  additional 
equation  is  required  to  model  the  dissipation  of  turbulence  e,  leading  to  the  so-called  k-e  equations  as  (see 
Justesen  1988  or  Wilcox  2000  for  the  complete  derivation): 

dt  dz 
and 

ds  __  d 
dt  dz 

The  use  of  Cl  in  the  dissipation  terms  follows  Justesen  (1988)  and  is  a  simplified  method  of  reducing 
the  dissipation  rather  than  using  empirical  functions  (see  Hwang  and  Jaw  (1998)  for  a  description  of  these 
various  functions). 

The  constants  used  in  Eqs.  (10),  (13),  and  (14)  are  often  found  from  direct  numerical  simulations  with 
the  typical  values  given  in  Table  1. 


i;+- 


'ej 


de 

dz 


(14) 


r 


o, 


\ 


y'^<Xk)dz\ 


dk 


+  u, 


_ 

.  ^  J 


(13) 


Table  1  —  Constants  Used  in  the  k-e  Model 


C; 

c, 

Cs 

0.09 

1.44 

1.92 

1.0 

1.3 
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Upon  determining  it  and  e,  the  turbulent  eddy  viscosity  is  determined  by 

y,=C,— .  (15) 

e 

It  is  known  that  the  ks  model  has  particular  difficulty  in  complex  flows  and  those  with  separation 
points  and  adverse  pressure  gradients  (Rodi  and  Scheuerer  1986,  Pope  2000,  and  Wilcox  2000).  Hence,  it 
should  be  used  in  direction-reversing  Iwundaiy  layer  flows  with  caution. 

Two-Equation  k-ca  Model 

Like  the  k-s  model,  the  k-(o  model,  which  is  known  to  be  more  robust  than  the  k-e  in  regions  of  adverse 
pressure  gradients  (e.g.,  Wilcox  2000),  is  a  two-equation  turbulence  closure  scheme  except  that  the 
second  equation  models  the  specific  dissipation  rate  <u  (=  elk),  rather  than  the  dissipation  rate  e. 

The  model  equations  (Wilcox  1988)  are 


and 


dt 


♦  \dk 
v-¥a  o,  \— 
‘Idz 


+  Ot 


r^f 
<  &  ) 


-J3*ka>, 


da> 

IT 


and  the  coefficients  are  given  in  Table  2. 


Table  2 — Constants  Used  in  the  k-o)  Model 


a 

* 

a 

B 

a 

0.5 

0.5 

3/40 

9/100 

5/9 

NUMERICAL  METHODS 
Finite  Difference  Equations 


(16) 


(17) 

(18) 


The  equations  for  the  boundary  layer  flow  and  the  various  eddy  viscosity  approaches  are  solved 
numerically  via  an  implicit  algorithm  using  2“*  order  spatial  differences  using  a  logarithmically  spaced 
grid.  The  implicit  Euler  method  is  given  as  (using  deflect  velocity  as  an  example) 


a,,  ,,n+l  ,,n 

dt  A/ 


(19) 


where  the  superscript  is  the  time  level.  Second  order  spatial  differences  are  given  as 
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dz  VA  hi-i  +  hilk 


hj_i  hj 

^d,i+\  ,  “rf,M 

hi-l 


and 


- - 


Vi+^  L  ^  h-\  J  VA’ 

where  i  is  the  grid  point  location  and 


(20) 


(21) 


hi  =2;+,  -2,-. 


(22) 


Systems  of  equations  are  developed  upon  substituting  Eqs.  (19)  through  (22)  for  the  various  variables 
into  the  model  equations. 

Boundary  Conditions 

In  order  to  folly  specify  the  problem,  boundary  conditions  are  needed  for  Eqs.  (10),  (13),  (14),  (16), 
and  (17).  For  foe  velocity,  foe  boundary  condition  imposed  at  foe  bed  is  no  slip  such  that  foe  deflect 
velocity  at  foe  roughness  level  (foe  first  grid  point)  is 


[30 


A 


j 


=  -Uit), 


(23) 


where  ku  is  foe  Nikuradse  equivalent  sand  roughness  and  zo  =  itj/JO  is  foe  roughness  length  where  foe 
fluid  velocity  goes  to  zero.  The  velocity  gradient  is  assumed  to  vanish  at  foe  top  of  foe  domain,  leaving 


du, 


dz 


^(D,r)  =  0. 


(24) 


The  boundary  conditions  for  k  are 


and 


dz 


(AO=0. 


(25) 


(26) 


The  lower  boundary  condition  on  k  stems  from  foe  fact  that  foe  turbulent  velocity  fluctuations  need  not  go 
to  zero  on  a  hydraulically  rough  bed.  Moreover,  equilibrium  is  assumed  to  exist  between  foe  production 
of  turbulent  kinetic  energy  and  dissipation,  yielding  Eq.  (27);  see  Justesen  (1988).  Rodi  (1980)  gives  foe 
boundary  conditions  for  the  dissipation  equation  as 


£ 


& 
I  30’ 


KZq  ’ 


(27) 


and 
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The  bottom  boundary  condition  for  the  specific  dissipation  rate  is  (Wilcox  2000) 


(28) 


(29) 


Solutions  are  sensitive  to  the  value  of  a  Dirichlet  fiee  stream  boundary  condition  so  we  apply  a 
Neumann  condition  at  the  surface  similar  to  the  conditions  for  k  and  e. 


•^(A0  =  0-  (30) 

Free-Stream  Forcing 

The  models  can  be  forced  by  either  a  time  series  of  free  stream  velocity  measurements  or  any  analytic 
wave  representation  with  or  without  a  superimposed  current.  Analytic  wave  representations  used  in  this 
study  are  sinusoidal  arid  sawtooth  waves.  Sinusoidal  forcing  is  given  by 


(7(0  =  C/;nSin 


(31) 


where  U„  is  the  maximum  free  stream  flow  amplitude  and  T  is  the  sinusoidal  wave  period.  Sawtooth 
waves  can  be  forced  using  the  form  described  by  Fredsoe  and  Deigaard  (1992)  as 


where  a,  is  a  parameter  controlling  the  skewness  of  the  wave  and  T  is  the  sawtooth  wave  period. 
Boundary  Layer  Estimation 

The  boundaiy  layer  height  is  calculated  in  one  of  two  ways  as 

g-2  ^  =Ps  (hereinafter  referred  to  as  BLl),  (33) 

where  uj+U  is  the  depth-dependent  velocity.  Ps  is  the  ratio  of  the  velocity  to  that  of  the  free  stream  and 
is  usually  set  to  0.95  or  0.99;  the  location  where  this  occurs  is  defined  as  Ae  boundary  layer  thickness  d. 
Difficulties  can  arise  with  Eq.  (33)  because  the  free-stream  velocity  goes  through  zero  when  forced  by 
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oscillatory  and  aperiodic  direction-reversing  flows.  Hence,  a  second  option  for  calculating  the  boundary 
layer  thickness  is  to  determine  when  the  magnitude  of  the  gradient  of  passes  a  tolerance  threshold  as 


S  =  z 


dz 


>  (hereinafter  referred  to  as  BL2). 


(34) 


Starting  at  the  surface  and  moving  in  the  downward  direction,  the  location  where  the  magnitude  of  the 
velocity  gradient  meets  or  exceeds  specifies  the  height  of  the  boundary  layer.  This  method  can 
produce  realistic  boundary  layer  heights  even  during  times  of  zero  or  near-zero  flow  above  the  boundary 
layer. 


SOLUTION  PROCEDURE 


The  models  (laminar,  linear,  linear-exponential,  parabolic,  k,  k-e,  and  k-co)  have  been  coded  in  Matlab 
using  the  implicit  numerical  methods  mentioned  above.  Unfortunately,  the  models  represent  coupled 
systems  of  equations  relying  on  an  unknown  (n+1)  time  step.  The  implicit  solution  procedure  for  the 
turbulence  closure  models  requires  iteration  on  the  variables  {ua  kr,  uj,  k,  e;  or  uj,  k,  (o)  to  obtain  a 
converged  solution  at  each  time  step.  The  eddy  viscosity  models  do  not  require  iteration. 

Figure  1  shows  a  flow  chart  of  the  solution  procedure.  After  defining  the  variables  in  the  workspace 
and  the  initial  conditions,  the  deflect  velocity  is  solved  for  and  the  vertical  derivative  of  the  deflect 
velocity  determined.  In  the  case  of  the  linear,  linear-exponential  or  parabolic  model,  tihe  fnction  velocity 
is  next  found  and  used  to  update  the  turbulent  eddy  viscosity  for  the  next  time  step.  For  the  one-equation 
model,  the  k  value  at  time  n  +  1  is  determined  and  used  to  update  the  turbulent  eddy  viscosity  using  the 
length  scale  throu^  Eq.  (12).  These  steps  are  repeated  for  this  time  step  until  the  error  on  the  turbulent 
eddy  viscosity  is  smaller  than  an  iteration  tolerance.  In  the  two-equation  models,  the  k  equation  is  solved, 
followed  by  £  or  ct)  equation,  and,  finally,  the  turbulent  eddy  viscosity  is  calculated  through  Eq.  (15)  or 
(18).  Similar  to  the  one-equation  model,  these  steps  are  iterated  upon  until  the  error  in  the  turbulent  eddy 
viscosity  is  smaller  than  an  iteration  tolerance.  The  variables  are  then  updated  for  the  next  time  step,  the 
boundaiy  layer  thickness  and  bed  shear  stress  are  calculated,  and  the  next  time  step  is  initiated. 


Fig.  1  — Flow  diagram  of  boundaiy  layer  model  solution  procedure 
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MODEL  RESULTS 

All  seven  models  were  tested  against  each  other  for  two  generic  cases.  (The  linear-exponential  model 
is  not  shown  in  any  of  the  figures  because  of  the  similarity  to  the  linear  and  parabolic  eddy  viscosity 
model.)  The  first  case  (CASEl)  is  a  5-s  sine  wave  with  an  amplitude  of  0.8  m  s"’.  The  second  case 
(CASE2)  is  a  5-s  sawtooth  wave  with  an  amplitude  of  0.8  m  s'*  and  a.,  -  0.25.  In  both  cases,  =  IxlO"^ 
m,  the  domain  height  is  0.1  m,  the  grid  spacing  was  specified  as  logarithmic  (with  300  points),  and  2000 
time  steps  were  completed  for  each  of  four  wave  periods  (only  the  last  was  retained  for  comparison). 

Figure  2  shows  velocity  fi-om  four  phases  (at  0°,  45°,  90°,  and  135°)  of  the  CASEl  comparison.  At  0°, 
Ae  flow  has  reversed  fiom  negative  to  positive,  and  the  phase  lead  of  the  boundaiy  layer  flow  can  be  seen 
in  all  cases.  The  laminar  simulation  predicts  a  very  small  boundaiy  layer  height  and  the  largest  near  bed 
velocity  at  over  0.25  m  s"'.  The  other  simulations  have  similar  (to  each  other)  near-bed  velocity  maxima 
of  about  0.18  m  s'*  (0.22  m  s'*  for  k-s  model)  but  vary  higher  in  the  water  column.  The  krs  and  k-co 
velocity  predictions  decrease  with  elevation  more  rapidly  than  the  k,  linear,  or  parabolic  simulations  (note 
that  the  linear  and  parabolic  simulations  are  nearly  identical).  Also,  these  two  simulations  still  predict 
negative  flow  fi-om  about  0.02  m  to  0.045  m,  whereas  the  linear,  parabolic,  and  jt  simulations  do  not. 


Fig-  2  —  CASEl  velocity  profiles  for  various  wave  jAases:  laminar  (- . );  linear 

parabolic  (■  ■  ■  ■);  k  ( - );  k-e  (  ■  );  and  k-m  (mh)- 

At  45°,  when  the  flow  is  accelerating,  all  simulations  except  the  laminar  one  have  essentially  the  same 
structure  with  small  variations.  For  instance,  the  k,  linear,  and  parabolic  models  do  not  have  an  inflection 
point  near  z  =  0.02  m.  At  90°,  the  acceleration  is  0  and  the  free-stream  velocity  is  at  its  maximum.  Again, 
the  simulations  look  similar  with  only  subtle  differences.  For  instance,  the  k-e  and  k-co  velocities  increase 
with  respect  to  elevation  faster  than  do  the  other  simulations  (except  laminar).  Finally,  by  135°,  the  flow 
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is  decelerating.  Note  the  difference  in  vertical  structure  between  this  plot  and  that  for  45°.  Both  have  the 
same  fiee-stream  velocity,  but  the  flow  at  45°  is  accelerating  while  the  flow  at  135°  is  decelerating. 

Figure  3  shows  the  CASE2  profiles.  At  0°,  the  shapes  of  the  vertical  profiles  do  not  vary  much  from 
their  CASEl  counteiparts.  The  laminar  model,  as  before,  predicts  a  much  thinner  boundaiy  layer  and 
larger  near-bed  velocity.  At  45°  and  90°,  the  shapes  of  the  profiles  are  similar  to  those  of  CASEl ,  but  the 
velocity  magnitude  is  reduced  by  roughly  one  half  and  one  fourth,  respectively,  due  to  the  skewed  forcing 
of  the  sawtooth  wave.  The  wave  skewness,  however,  causes  the  velocity  maxima  at  135°  to  exceed  those 
in  CASEl  by  over  0.2  m  s'’. 


fig-  3  —  CASE2  velocity  profiles  for  various  wave  phases:  laminar  (-. . .);  linfar 

parabolic  (■■■■);  A:  ( - );  k-e  (  );  and  k-co  (m). 


Figure  4  shows  the  boundary  layer  estimations  based  on  the  two  methods  (BLl  and  BL2)  for  CASEl . 
The  advantage  of  BL2  (Fig.  4A)  over  BLl  (Fig.  4B)  is  clearly  seen  at  flow  reversals  when  the  free-stream 
velocity  passes  through  zero,  and  the  BLl  estimate  reaches  the  domain  height  While  the  boundary  layer 
estimation  using  the  velocity  gradient  appears  to  perform  better,  there  are  still  considerable  differences  in 
the  sh^  between  the  various  simulations.  Decreases  in  the  estimated  boundary  layer  height  are  seen 
for  the  and  k-co  simulations  near-flow  reversal,  while  decreases  are  seen  for  the  other  nonlaminar 
simulations  closer  to  peak  free-stream  velocity  magnitudes.  Other  estimations  for  the  BL2  method  can  be 
obtained  by  varying  the  cutoff  from  the  value  of  1  used  here. 


nl _ — T - 1 - - — I - - - - j _ I - 1 - j - r - 1 

15  15.5  16  16.5  17  17.5  18  18.5  19  19.5  20 

Time  (s) 


Fig.  4  —  CASEl  boundaiy  layer  estimations  for  BL2  (A)  and  BLl(B):  laminar  (• . );  linear 

(.. . — );  parabolic  . );  k  (  k-s  (  ■■);  k-m  (imv). 


Bed  shear  stress  estimates  for  CASEl  show  variation  over  the  period  and  between  the  different  model 
formulations  (Fig.  5).  For  instance,  all  stress  estimates  have  the  same  shape,  but  vaiy  in  magnitude  from 
a  maximum  of  nearly  1  kg  m’’  s‘^  for  the  laminar  simulation  (Fig.  5;  thin  dashed  curve)  to  just  over  4.0  kg 
m’’  s'^  for  the  linear  and  eddy  viscosity  model  simulations  (Fig.  5;  bold  dashed  curves).  These 
simulations  also  show  the  phase  lead  of  the  bed  shear  stress  with  respect  to  the  free-stream  velocity.  This 
is  apparent  because  the  free-stream  velocity  maxima  occur  at  16.25  and  18.75  s  while  the  two  peaks  in  the 
bed  shear  stress  occur  at  16.0  and  18.5  s  for  the  linear  and  parabolic  simulations  and  15.625  and  18.125 
for  the  laminar  simulation.  This  represents  a  phase  lead  of  0.25  s  (18°)  for  the  linear  and  parabolic 
simulations  and  0.625  s  (45°)  for  the  laminar  simulation.  The  45°  lead  for  the  laminar  simulation  exactly 
matches  the  phase  lead  predicted  by  boundary  layer  theory  (Eq.  (3))  for  sinusoidal  wave  forcing  (Stokes 
1845). 
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Fig.  5  —  CASEl  bed  shear  stress  estimate,  laminar  (• . );  linear  (.. . );  parabolic  (■  ■  ■  a); 


Bed  shear  stress  estimates  for  the  turbulent  closure  schemes  (Fig.  5;  solid  curves)  have  the  same  shape 
as  the  eddy  viscosity  estimates  but,  interestingly,  are  generally  smaller  in  magnitude.  The  phase  lead  of 
the  bed  shear  stress  of  15.8°,  12.2°,  and  18°  for  the  k,  k-e,  and  k-co  simulations,  respectively,  is  similar  to 
that  of  the  linear  and  parabolic  simulations  and  comparable  to  those  given  by  Fredsoe  and  Deigaard 

(1992,  pg.  32)  for  amplitude  Reynolds  numbers  ( RE  =  •^,  a  =  )  used  here.  The  inflection  in 

the  k-e  model  bed  shear-stress  estimate  is  likely  attributed  to  the  flow  reversal. 

MODEL-DATA  COMPARISON 

Jensen  et  al.  (1989)  collected  depth-dependent  flow  velocity  measurements  in  a  U-shaped  oscillatory 
water  tunnel  for  smooth  and  rough  bed  cases  using  a  laser  Doppler  anemometer  (LDA).  The  tunnel  was 
10  m  long,  0.29  m  wide,  and  0.3  m  deep  (0.28  m  in  smooth  wall  cases).  The  models  will  be  tested  against 
the  J^en  et  al.  (1989)  TestlO  and  Testl3.  Testl3  was  forced  with  a  sinusoid  of  9.72  s,  a  maximum 
velocity  of  2  m  s  and  a  roughness  of  Kn  =  0.84  mm.  Testl  0  was  forced  with  the  same  conditions  except 
the  bed  was  made  smooth  by  the  installation  of  PVC  plates.  Because  the  models  require  a  zq  (=  Kfi/30) 
value,  a  smooth  bed  was  approximated  by  using  =  1  xlQ-^  mm  (four  orders  of  magnitude  smaller  than 
the  rough  bed  case  and  one  order  of  magnitude  smaller  than  the  largest  clay  particle  size  on  the 
Wentworth  scale).  Similar  to  the  generic  cases,  2000  time  steps  per  wave  period  were  used  and  the 
number  of  grid  points  was  300. 


Wave  Bottom  Boundary  Layer  Models  for  Smooth  and  Rough  Beds  _  13 

Figures  6  and  7  show  the  velocity  comparison  for  eveiy  15  degrees  of  phase  for  the  first  half  of  the 
wave  for  the  simulations  (solid  lines)  and  the  Jensen  et  al.  (1989)  Testl3  (dots).  Each  successive  profile 
is  ofifeet  by  2  m  s'*  (denoted  by  the  zeroes  in  the  plot).  The  panels  in  Fig.  6  are  A)  laminar,  B)  linear,  and 
C)  parabolic.  Similarly,  the  panels  in  Fig.  7  are  A)  k,  B)  k-e,  and  C)  k-o>.  The  laminar  simulation  has 
difficulty  matching  the  observed  profile  shape  at  all  phases  of  wave  motion; 


Fig.  6  —  Rough  wall  (Testl3),  experimental  data  (.),  model  (solid  curves).  A)  laminar,  B)  Imear,  and 
C)  parabolic  eddy  viscosity. 
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Fig.  7—  Rough  wall  (Testl3),  experimental  data  (.),  and  model  (solid  curves):  A)  k,  B)  k-E ,  and  C) 
k-w  turbulence  models. 


smulation  generally  has  the  correct  velocity  profile  shape  and  magnitude  in  the  upper  portion 
of  fte  domain.  Near  the  bed,  however,  the  k-e  simulation  only  matches  the  shape  of  the  observations 
underpredicting  the  velocity  magnitude  at  nearly  all  phases  of  the  wave.  The  rest  of  the  simulations  are 
very  similar  and  qualitatively  match  the  data  (quantitative  comparisons  are  given  below). 

Figures  8  and  9  show  the  velocity  comparisons  to  the  data  of  Jensen  et  al.  (1989)  TestlO  (layout  is  the 
same  as  m  Figs.  6  and  7).  The  results  are  similar  to  the  rough  bed  test  in  that  the  laminar  simulation  is  the 
worrt  performer  and  the  k-e  simulation  underpredicts  the  near-bed  velocity  at  nearly  all  phases  of  wave 
motion.  Figures  6  through  9  show  that  the  simulations  are  qualitatively  similar  regardless  of  whether  or 
not  they  are  applied  to  rough  or  smooth  beds. 


Elevation  (m) 
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Fig.  9 —  Smooth  wall  (Test  10),  experimental  data  (.),  and  model  (solid  curves):  A)  k,  B)  k-e ,  and 
C)  k-<o  turbulence  models. 
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A  quantitative  measure  of  the  predictive  capability  of  the  models  can  be  determined  fix>m  the 
normalized  time  dependent  error 


_  ~“l 


^  data 


mode/  J 

(0 


(35) 


where  is  the  time-dependent  vertical  variance  in  data.  values  were  calculated  using  the  model 

grid  points  closest  to  the  measurement  locations.  This  distance  was  typically  on  the  order  of  1x10'’  m  or 
less  except  in  die  upper  few  centimeters  of  the  domain  where  the  maximum  distance  was  1.6x10'^  m.  The 
smooth-  and  rou^-bed  comparisons  show  that  the  linear  model  was  most  accurate  (although  the 
parabolic,  k,  and  k-co  models  had  similar  errors)  while  the  laminar  model  was  least  accurate  (see  Fig.  10; 
laminar  errors  extend  beyond  the  axis  range  shown).  The  k-e  model  had  normalized  errors  that  were  up  to 
six  times  as  large  as  those  of  the  k  and  k-co  model.  The  largest  errors  are  observed  between  0®  and  20° 
and  between  160°  and  180°  near  flow  reversals.  Away  from  flow  reversals  (adverse  pressure  gradients), 
all  of  the  models  except  the  laminar  and  k-e  perform  nearly  equally  for  both  the  rough  bed  (Fig.  10  A)  and 
smooth  bed  (Fig.  lOB)  tests.  In  general,  the  best  performers  for  both  tests  were  the  linear  and  parabolic 
eddy  viscosity  simulations  and  the  k-co  model.  Note  the  difference  in  scales  on  Figs.  lOA  and  lOB, 
implying  the  models  as  a  whole  perform  more  poorly  (by  up  to  a  factor  of  3)  for  the  smooth-bed 
simulation. 


Fig.  10  —  Normalized  errors  for  rough-wall  Testl3  (A)  and  smooth-wal!  TestlO  (B):  laminar  (gray  *-); 
linear  {♦);  parabolic  (square);  k  (diamond);  k-e  (o);  and  k-co  (x). 
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DISCUSSION  AND  CONCLUSIONS 

The  performance  of  seven  boundaiy  layer  models  has  been  investigated  for  two  generic  cases  (sine 
wave  and  sawtooth  wave)  and  data  collected  in  a  laboratory  U-shaped  oscillatory  water  tunnel.  The  eddy 
viscosity  boundary  layer  models,  where  the  vertical  profile  of  the  turbulent  eddy  viscosity  is  assumed  to 
have  a  specific  shape  were  similar  in  their  predictive  capability.  Their  similarity  stems  from  the  fact  that 
within  the  boundary  layer,  the  linear,  linear-exponential,  and  parabolic  forms  for  the  eddy  viscosity  are 
also  similar.  The  shapes  of  the  ed<fy  viscosity  profile  for  these  varies  outside  the  boundaiy  layer,  but 
because  there  is  little  shear  in  the  water  column  at  these  elevations,  the  differences  have  little  effect  on  the 
predicted  velocity  profiles.  The  effect  of  a  superimposed  current  is  not  included  in  the  eddy  viscosity 
models  used  here,  but  can  be  included  through  a  modified  eddy  viscosity  by  either  including  the  current 
effect  in  the  bed  shear  stress  estimate  (e.g..  Grant  and  Madsen  1979)  or  in  a  region-dependent  (i.e.,  in  the 
boundary  layer  or  above)  eddy  viscosity  model  (e.g.,  Fredsoe  and  Deig^d  1992). 

The  accuracy  for  the  turbulence  closure  schemes  cannot  be  determined  from  the  generic  cases;  rather, 
die  attempt  was  to  see  how  these  models  compared  to  the  simplified,  less  numerically  intensive  eddy 
viscosity  models.  Based  on  these  generic  cases,  it  does  not  appear  that  the  turbulence  closure  models 
yield  velocity  profiles  that  differ  much  from  the  eddy  viscosity  models.  There  was  discrepancy,  however, 
between  the  predicted  bed  shear  stress  between  the  two  model  types,  with  the  eddy  viscosity  models 
predicting  larger  values  than  the  turbulence  closure  schemes.  These  differences  stem  from  the  fact  tfiat 
the  eddy  viscosity  models  have  a  predefined  shape  for  the  eddy  viscosity  whereas  the  turbulence  closure 
models  do  not. 

The  comparison  to  laboratory  data  showed  how  accurately  the  eddy  viscosity  and  turbulence  closure 
models,  with  the  exception  of  the  k~e  model,  predicted  the  observations.  For  instance,  the  k  and  k-<o 
models  performed  well  at  all  phases  of  the  wave.  The  eddy  viscosity  models  also  compared  well  to  the 
laboratory  data  for  smooth  and  rough  beds  except,  at  times,  in  regions  of  flow  reversal.  The  consistently 
larger  errors  associated  with  the  k-e  model  result  from  the  inability  of  the  model  to  compensate  at  later 
phases  for  the  overprediction  during  flow  reversal.  It  has  been  suggested  and  shown  previously  (Rodi  and 
Scheuerer  1986  and  Wilcox  1988)  that  the  k-e  model  can  peiTorm  poorly  during  flow  reversal  (in  regions 
with  an  adverse  pressure  gradient)  as  was  observed.  A  typical  workaround  that  was  shown  here,  use  of 
the  k-Q}  model,  has  better  predictive  capability  during  these  times  (Wilcox  1 988). 

The  WBBL  models  presented  here  are  a  means  by  which  vertical  velocity  profiles  and  bed  shear 
stresses  cm  be  generated  for  a  variety  of  analytic  wave  forms  and  free-stream  velocity  time  series.  The 
purpose  of  the  models  is  to  provide  a  simple  means  to  estimate  these  parameters  under  a  variety  of 
conditions  with  the  hopes  of  improving  the  knowledge  of  the  fluid  flow  within  the  boundaiy  layer  and 
potentially  enhancing  sediment  transport  predictive  capability.  Based  on  this  analysis,  it  does  not  appear 
that  the  extra  computation  time  for  the  turbulence  closure  schemes  (six  times  slower  than  the  eddy 
viscosity  models  for  the  number  of  time  steps  used  here)  affords  improvement  in  predictive  capability.  In 
conditions  with  much  higher  Reynolds  numbers  (i.e.,  cases  with  much  larger  velocities),  stronger  adverse 
pressure  gradients,  and/or  2D  or  3D  simulations,  the  use  of  the  turbulence  closure  models  may  be  more 
appropriate. 
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Appendix 
MATLAB  CODES 


The  Matlab  (written  using  Matlab  6.0)  codes  shown  in  Table  A1  are  included  on  the  accompanying 
disk  for  convenience  to  the  reader.  The  codes  are  heavily  commented  to  explain  the  procedure 
throughout.  The  included  setvars.m  file  will  perform  the  simulation  referred  to  as  CASEl  for  linear  eddy 
viscosity.  Every  variable  is  described  after  it  is  defined  and  the  recommended  values  are  used.  As  an 
example,  at  the  Matlab  prompt,  type:  setvars(‘caser).  Then  type  evmodel(‘caser).  To  see  a  quick  plot 
of  the  velocity,  load  easel  .mat  and  type  plotmodel. 


Table  A1  —  Matlab  Codes  Used  for  Modeling  the  WBBL 


setvars.m 

Sets  up  the  variables  and  constants  for  a  model  run.  Used  for  all  4  modelm 
files 

evmodel.m 

The  WBBL  model  for  the  laminar,  linear,  linear-exponential,  and  parabolic 
eddy  viscosity  model  cases 

kmodel.m 

The  WBBL  model  using  the  one-equation  turbulence  model 

kemodeLm 

The  WBBL  model  using  the  two-equation  k-s  turbulence  model 

kwmodeLm 

The  WBBL  model  using  the  two-equation  k-co  turbulence  model 

plotmodeLm 

Plots  the  velocity  profiles  for  the  last  wave  in  the  simulation  at  45-deg  intervals. 

